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O ; ABSTRACT 

(N 

p ' The gravitational wave signal generated by global, nonaxisymmetric shear flows 

■^ ' in a neutron star is calculated numerically by integrating the incompressible 

00 ■ Navier-Stokes equation in a spherical, differentially rotating shell. At Reynolds 

numbers i?e ^ 3 x 10^, the laminar Stokes flow is unstable and helical, oscillating 
Taylor-Gortler vortices develop. The gravitational wave strain generated by the 
fT^ ■ resulting kinetic-energy fluctuations is computed in both + and x polarizations 

^ ■ as a function of time. It is found that the signal-to-noise ratio for a coherent, 10*- 

'sj- ' s integration with LIGO II scales as 6.5(r2^,/10^ rad s~^)^' ^ for a star at 1 kpc with 

^ ■ angular velocity il^. This should be regarded as a lower limit: it excludes pressure 

O ■ fluctuations, herringbone flows, Stuart vortices, and fully developed turbulence 

O : (for Re > 10^). 
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1. Introduction 

Gravitational radiation from linear, global, fluid oscillations in neutron stars has been 
studied extensively, e.g. radiation-reaction-driven r-modes (Andersson 1998; Levin & Ushomirsky 
2001), and two-stream superfluid oscillations with entrainment (Andersson 2003; Prix et al. 
2004). However, gravitational radiation is also emitted by nonlinear, global, fluid oscilla- 
tions. Specifically, it is well known that a viscous Navier-Stokes fluid inside a spherical, 
differentially rotating shell — spherical Couette flow (SCF) — undergoes sudden transitions 
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between states with different vortex topologies, which are generally nonaxisymmetric, due to 
shear instabilities (Marcus & Tuckerman 1987a,b; Junk & Egbers 2000). Recent simulations 
of superfiuid SCF, based on the Hall-Vinen-Bekarevich-Khalatnikov (HVBK) theory of a 
^5*0 neutron superfiuid, confirm that such transitions to nonaxisymmetric flows also occur in 
neutron stars (Peralta et al. 2005, 2006). 

In this Letter, we compute the gravitational wave signal generated by nonaxisymmetric 
flows in young, rapidly rotating neutron stars. There are at least two astrophysical scenarios 
in which these develop. First, if the crust and core of the star are loosely coupled, as indi- 
cated by observations of pulsar glitches (Shemar & Lyne 1996; Lyne et al. 2000), differential 
rotation builds up as the crust spins down electromagnetically, until the core undergoes a 
transition to the Taylor-Gortler vortex state, a helical, oscillating flow that is one of the 
last stages in the laminar-turbulent transition in SCF as the Reynolds number increases 
(Nakabayashi 1983; Nakabayashi & Tsuchida 1988). Second, if the crust precesses, while the 
core does not, the rotation axes of the crust and core misalign, inducing a rich variety of 
flow patterns and instabilities, like time-dependent shear waves, leading to fully developed 
turbulence (Wilde & Vanyo 1995). We treat the former scenario here and postpone the latter 
to a separate paper. 

In §2, an SCF model of a differentially rotating neutron star is proposed. In §3, we ex- 
plore the transitional dynamics and topology of the flow. In §4, we calculate the gravitational 
wave signal and its frequency spectrum. 



2. Global SCF model and numerical method 

We consider an idealized, two-component model of a neutron star, in which a solid crust 
is loosely coupled to a differentially rotating fluid core. Specifically, we consider the motion 
of a neutron fluid within a spherical, differentially rotating shell in the outer core of the star, 
where the density p lies in the range 0.6 < p/ Pc < 1-5, Pc = 2.6 x lO^^gcm"^ (Sedrakian & 
Sedrakian 1995). The inner (radius -Ri) and outer (radius R2) boundaries rotate at angular 
frequencies Qi and ^2, respectively, about the z axis. As the star is strongly stratified, the 
shell is thin, with R2 - Ri < O.lRi (Abney & Epstein 1996; Levin & D'Angelo 2004). The 
rotational shear is sustained by the vacuum-dipole spin-down torque, violent events at birth 
(Dimmelmeier et al. 2002; Ott et al. 2004), accretion in a binary (Fujimoto 1993), neutron 
star mergers (Shibata & Uryii 2000), or internal oscillations like r-modes (Rezzolla et al. 
2000; Levin & Ushomirsky 2001). 

We describe the fluid by the isothermal Navier-Stokes equation, which, in the inertial 
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frame of an external observer, takes the form 

|^ + (vV)v = -^ + z/VV + V<l> (1) 

in the incompressible limit V • v = 0, where v is the fluid velocity, v is the kinematic 
viscosity, p is the density, p is the pressure, and $ is the Newtonian gravitational potential. 
Henceforth, $ is absorbed into p by replacing p — p$ with p, such that pressure balances 
gravity in the stationary equilibrium (v = 0) and there are no gravitational forces driving the 
flow in the incompressible limit when the stationary equilibrium is perturbed. The Reynolds 
number and dimensionless gap width are defined by Re = QiRf/u and 6 = {R2 — Ri)/Ri] 
the rotational shear is AQ = ^^2 — fii. Equation (1) is solved subject to no-slip boundary 
conditions. 

Superfluidity plays a central role in the thermal (Yakovlev et al. 1999) and hydrody- 
namic (Alpar 1978; Reisenegger 1993; Peralta et al. 2005) behavior of neutron star interiors. 
However, we assume a viscous fluid here. Counterintuitively, this is a good approximation 
because Re (~ 10^^) is very high (Mastrano & Melatos 2005). It is known from terrestrial 
experiments on superfiuid '^He that, at high Re, quantized vortices tend to lock the superfluid 
to turbulent eddies in the normal fluid by mutual friction (Barenghi et al. 1997), so that the 
superfluid resembles classical Navier-Stokes turbulence. Moreover, in HVBK simulations of 
SCF, it is observed that the global circulation pattern of the superfluid does not differ much 
from a Navier-Stokes fluid at Re ^ 250 (Henderson & Barenghi 2004; Peralta et al. 2005). 
Finally, a viscous interior is expected in newly born neutron stars, which are hotter than the 
superfluidity transition temperature (Andersson et al. 1999). 

By assuming a viscous fluid, we omit from (1) the coupling between quantized vortices 
and the normal fluid (Hall & Vinen 1956), entrainment of superfluid neutrons by protons 
(Comer 2002), and pinning of quantized vortices in the inner crust, at radius Ri [Baym et al. 
(1992); cf. Jones (1998)]. We also neglect vortex pinning in the outer core for simplicity, 
even though evidence exists that it may be important when magnetic flelds are included. 
Speciflcally, three-fluid models of post-glitch relaxation based on the core dynamics favor 
vortex pinning at the phase separation boundary of the core (Sedrakian & Sedrakian 1995), 
e.g. due to the interaction between vortex clusters and the Meissner supercurrent set up by 
the crustal magnetic fleld at the phase boundary (Sedrakian & Cordes 1999). Other models 
analyze interpinning of proton and neutron vortices in the core (Ruderman 1991) and its 
effect on precession (Link 2003). In the context of the hydrodynamic model in this paper, 
pinning effectively increases the viscosity of the core, reducing Re. 

We use a pseudospectral collocation method to solve equation (1) (Bagchi & Balachan- 
dar 2002; Giacobello 2005). The equations are spatially discretized in spherical polar coor- 
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dinates {r,6,(f)), using restricted Fourier expansions in 6 and and a Chebyshev expansion 
in r. The solution is advanced using a two-step fractional-step method, which is accurate 
to second order in the time-step (Canuto et al. 1988). We limit ourselves to narrow gaps 
6 < 0.1, which are computationally less expensive and for which more experimental stud- 
ies are available for comparison (Yavorskaya et al. 1977; Nakabayashi 1983). Recently, we 
performed the first stable simulations of superfluid SCF using this method (Peralta et al. 
2005). 



3. Nonaxisymmetric spherical Couette flow 

The bifurcations and instabilities leading to transitions between SCF states are con- 
trolled by three parameters: S, Re, and AQ. Additionally, the history of the flow influences 
its post-instability evolution and the final transition to turbulence (Wulf et al. 1999; Junk & 
Egbers 2000). For example, in experiments on Taylor-Gortler vortices (TGV), the final state 
depends on dQi/dt, i.e. on the time to reach the critical Re relative to the viscous diffusion 
time (Nakabayashi & Tsuchida 1995). In general, four SCF states can be distinguished: (i) 
a laminar basic flow, (ii) a toroidal or helical Taylor-Gortler flow, (iii) a transitional flow 
with nonaxisymmetric, oscillating TGV, and (iv) fully developed turbulence (Yavorskaya 
et al. 1975; Nakabayashi 1983). In a neutron star, where Re ~ 10^^, the state is probably 
turbulent. However, we are restricted by our computational resources to simulate regimes 
(i) and (ii), with 1 x 10^ < i?e < 3 x 10^. 



3.1. TGV transition: initial conditions 

Experimentally, for 6 = 0.06, TGV are obtained by quasistatically increasing Qi from 
zero to give Re ~ 3300 (Nakabayashi & Tsuchida 2005). This procedure can be painfully 
slow, because a steady state must be reached at each intermediate step. Numerically, we 
circumvent it by following Li (2004) and introducing a nonaxisymmetric perturbation with 
azimuthal wave number rria of the form 

4ei{r - Ri){r - R2) ,„ ,. ,_. 

■" = ^ cosx{0,(j)), (2) 



V0 = — sm 




5^Rl 
TTJr - Ri){r - R2 
2(i?2 - Riy 



sinx(^,0), (3) 



with X = 7r[l — -R2(7r/2 — 6)/{R2 — Ri) — 0.4 sm{ma(f))]. This is done as follows. Before intro- 
ducing the perturbation, a steady state for Re = 2667 is obtained (just below the critical Re 
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where TGV emerge experimentally). Then, the Reynolds number is raised instantaneously 
to Re = 3300, and we continue by adding (2) and (3), with amplitude ci ~ 10~® QiRi, to the 
numerical solution at each time-step until a viscous diffusion time td = {R2 — RiY/u elapses. 
We then stop adding the perturbation and the ffow is left to evolve according to (1), until a 
final steady state is reached. More than one perturbation can lead to the same final state; 
some authors add Gaussian noise (Zikanov 1996), although this affords less control over the 
wavenumber and the number of vortices excited. 

The perturbation excites TGV by shedding vorticity from the inner sphere (Li 2004). 
Perturbations with 2 < ma < 5 are explored, supplementing experiments with nia = 3 
(Nakabayashi 1983). For numerical simplicity, we limit ourselves to the case where only the 
inner sphere rotates. TGV with both spheres rotating are equally possible and have been 
observed experimentally; they exhibit additional twisting (and increased nonaxisymmetry) 
in the helical vortices (Nakabayashi & Tsuchida 2005), so the gravitational wave strain we 
compute in §4 is a lower limit. 



3.2. Flow topology 

Figure la shows a kinetic-energy-density isosurface (0.11pi?ff2^) for the fully developed 
TGV state rria = 3, 6 = 0.06, Re = 3300. Its nonaxisymmetry is apparent in the striated 
equatorial bands, whose inclination with respect to the equator varies with longitude from 
zero to ~ 3 deg (Nakabayashi 1983). 

A popular way to classify complex, three-dimensional flows is to construct scalar in- 
variants from Aij = dvi/dxj (Chong et al. 1990; Jeong & Hussain 1995). Specifically, the 
discriminant Da = Qa + '27Ra/^i with Ra = -det{Aij) and Qa = (A^ - AijAji)/2, distin- 
guishes between regions which are focal {Da > 0) and strain-dominated {Da < 0). Figure 
lb plots the Da = 10~^ Qf isosurface in color: yellow regions are stable focus/stretching, 
^ blue regions are unstable focus/contracting. The filaments coincide with the TGV. One 
circumferential vortex and rria = 3 helical vortices lie in each hemisphere. We also plot the 
isosurface Da = — 10^^f2p red and green regions have stable node/saddle/saddle and un- 
stable node/saddle/saddle topologies respectively (Chong et al. 1990). The helical vortices 
span ~ 40 deg of latitude and travel in the direction, with phase speed ~ OA8Q1R1. 

The TGV state is oscillatory, which is important for the gravitational wave spectrum. 



^Trajectories are repelled away from a fixed point along a real eigenvector of Aij (stretch) and describe 
an inward spiral when projected onto the plane normal to the eigenvector (stable focus). 
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Its periodicity is evident in Figure Ic, where vg at an equatorial point is plotted versus time 
for 2 < nia < 5. For rria = 3, vg oscillates with period 4.4^2^"^, the time for successive helical 
vortices to pass by a stationary observer, in accord with experiments (Nakabayashi 1983). 
Instantaneous streamlines are drawn in three meridional planes in Figure Id, highlighting the 
nonaxisymmetry. One obtains a similar TGV state for 6 = 0.14, = 0.18 (Sha & Nakabayashi 
2001; Li 2004). 



4. Gravitational wave signal 

The metric perturbation hj^ = {AG/rc^) J rf^x Tij in the transverse-traceless gauge can 
be calculated using Einstein's quadrupole formula (Misner et al. 1973) 

|2 



h¥ = —r^^ / c^'x |v|2 {x,Xj - 5i}^ ) , (4) 



where the integral is over the source volume, r is the distance to the source, and Tij = 
P'y'^ViVj+p6ij is the stress-energy of a Newtonian fluid. Note that we approximate Too = P7^ ~ 
p(l + |vp/c^) in (4), and we omit a thermal-energy contribution oc 6p (~ p|vp, by Bernoulli's 
theorem), which we cannot calculate with our incompressible solver. Far from the source, 
we write hj,^ = h+ef, + hy^ef,, with the polarizations deflned by e,t, = — ej; = e^,, = e^„ = 1 
for an observer on the x axis. 

In Figure 2a, we plot the h^ and /ix polarizations versus time, as seen by an observer 
on the X axis, for 2 < irta < 5 (the signal is ~ 3 times weaker on the z axis, where the 
north-south asymmetry is not seen in projection). The amplitude is greatest for rria = 5 
but depends weakly on m^. Importantly, /i+ and h^ are 7i/2 out of phase for all ma, and 
we find |/i+| <^ \hx\ for even nia- These two signatures offer a promising target for future 
observations. The period (~ 0.8i7f ^) is similar for 3 < ma < 5, where a maximum of three 
helical vortices are excited (Li 2004) (we flnd that the torque at R2 oscillates with the same 
period). This too is good for detection, because several modes are likely to be excited in a 
real star. The period arises because the isosurface in Figure la forms a pattern with six-fold 
symmetry when projected onto the y-z plane, with fundamental period 4.4^2^^^ 

In Figure 2b, we present the frequency spectra of the + and x polarizations. The two 
peaks, at / ~ 1.2i7i and 2.4r2i, have full-width-half-maxima of ~ 300 Hz. This is caused by 
the subharmonics evident in Figure Ic, which arise because the northern and southern helical 
vortices start at unequal and variable longitudes. In addition, the peaks for even and odd 
nia CL'>^G displaced by ~ 500 Hz, and the primary peak for nia = 4 zs split. These two spectral 
signatures are a promising target for future observations. The even-odd displacement may 
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arise because the phase speeds of the vortex patterns for rria = 2,3 differ by ^ 0.0ir2ii?i (Li 
2004). For large /, we find /i+,x(/) oc /"^Z^. 

The squared signal-to-noise ratio can be calculated from {S/Ny = 4: f^ df\h{f)\'^/Sh{f) 
(Creighton 2003), where we take Sh{f) = IQ-^^f/O.GkRz)^ (0.2 < / < 3 kHz) for a 10^-s 
integration with LIGO II, if the frequency and phase of the signal are known in advance 
(Brady et al. 1998). For a star with f]i/27r = 600 Hz, at r = 1 kpc, we find S/N = 0.21 
{rria = 2), 0.12 {rria = 3), 0.023 {rria = 4), and 0.22 {rria = 5), although there will be some 
leakage of signal when the coherent integration is performed, because the spectrum is not 
monochromatic. For m^ < 5, S/N increases (decreases) as rria odd (even) increases. We 
also find S/N oc Q^ , which implies that the fastest millisecond pulsars {fii/27i ~ 1 kHz) 
and newly born pulsars with Ar2/r2i ~ 0.1 are most likely to be detected. The flow states 
leading to turbulence at Re ^ 10^, e.g. shear waves, Stuart vortices, and "herringbone" 
waves (Nakabayashi & Tsuchida 1988), make contributions of similar order to the TGV 
signal in Figure 2. 

The Reynolds number in the outer core of a neutron star. Re = 10^''(p/10^^gcm^^)~^ 
(T/10^K)2(f]^/10'^rads"^) (Mastrano & Melatos 2005), where T is the temperature, typ- 
ically exceeds the maximum Re in our simulations (due to computational capacity). For 
Re ^ 10^, the flow is turbulent and therefore axisymmetric when time- averaged — but not 
instantaneously. In Kolmogorov turbulence, the characteristic flow speed in an eddy of 
wavenumber k scales as k~^^^, the turbulent kinetic energy scales as k~^^^, and the turn- 
over time scales as /c^^^^. From equation (4), the rms wave strain scales as /c~^/^ J dk k~^^^ 
after averaging over Gaussian fluctuations, and is dominated by the largest eddies (where 
most of the kinetic energy also resides). We therefore conclude that, even for Re ^ 10^, the 
largest eddies, which resemble organized structures like TGV, dominate h^ and h^, not the 
isotropic small eddies, but the signal is noisier. Hence, counterintuitively, we predict that 
hotter neutron stars (with lower Re) have narrower gravitational wave spectra than cooler 
neutron stars, ceteris paribus. 

We ignore compressibility when computing Tqo {Sp = 0) and solving (1) numerically 
by pressure projection. Yet, realistically, neutron stars are strongly stratifled (Abney & 
Epstein 1996), conflning the meridional flow in Figure Id into narrow layers near r = Ri, 
i?2- Stratiflcation can be implemented crudely by using a low-pass spectral fllter to artiflcially 
suppres Vr (Don 1994; Peralta et al. 2006). We postpone this to future work. 

We acknowledge the computer time supplied by the Australian Partnership for Advanced 
Computation (APAC). We thank Prof. Li Yuan, from the Chinese Academy of Sciences, for 
very helpful discussions. 
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Fig. 1. — Taylor-Gortler vortices for nia = 3, 5 = 0.06, at t = lOOf^jf^. (a) Isosurface 
of kinetic energy: p|vp = O.llpi^ffi^ (b) Isosurfaces of velocity-gradient discriminant: 
Da = 10"^ fi^, stable focus/stretching (yellow) and unstable focus/contracting (blue); Da = 
— 10^^i7^, stable node/saddle/saddle (red) and unstable node/saddle/saddle (green), (c) 
Time history of f g at r = 1.06, 6 = 7r/2, = for rria = 2 (red), 3 (blue), 4 (green), and 5 
(black), (d) Streamlines obtained by integrating the in-plane components of v in the planes 
(j) = (left), (f) = 7r/2 (middle), and = 27r/9 (right). 
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Fig. 2. — Gravitational waves from Taylor-Gortler vortices in a neutron star with Qi/27i = 
600 Hz and r = 1 kpc. (a) Signal h^^y.(t), and (b) frequency spectrum /i+ xl/), for rria = 2 
(red), 3 (blue), 4 (green), and 5 (black). Solid and dashed curves correspond to the + and 
X polarizations respectively, as measured by an observer on the x axis. 



